#ifndef __PSI_SW_BLOCK_3_KERNEL_H_
#define __PSI_SW_BLOCK_3_KERNEL_H_

namespace new_hope_3 {

const int constant_query_max = 60000;
const int alphabet_len = 25;

/* using regular includes didnt work - dont know why */
#include "include_current"


/** Smith-Waterman gap open cost */
__constant__ int gap_open;

/** Smith-Waterman gap extend cost */
__constant__ int gap_extend;

/** Ordinary score matrix expanded to profile matrix. */
__device__ int profile_matrix[constant_query_max][alphabet_len];


__device__ void single(int& h_max, int& f_top, int& h_top, int& e, int& h, int& h_topleft, const int& score){

	int f = max(f_top - gap_extend, h_top - gap_open);
	e = max(e - gap_extend, h - gap_open);
	h = max(0, max(e, max(f, h_topleft + score)));

	f_top = f;
	h_topleft = h_top;
	h_top = h;

	h_max = max(h, h_max);
};

__device__ void do_8step(short2* eh_tmp, int& h_topleft, int& h_top_0, int& f_top_0, int& h_top_1, int& f_top_1, int& h_top_2, int& f_top_2, int& h_top_3, int& f_top_3, int& h_top_4, int& f_top_4, int& h_top_5, int& f_top_5, int& h_top_6, int& f_top_6, int& h_top_7, int& f_top_7, int& h_max, const int& score_0, const int& score_1, const int& score_2, const int& score_3, const int& score_4, const int& score_5, const int& score_6, const int& score_7, const int& offset){
	short2 tmp = eh_tmp[offset];

	int e_left = tmp.x;
	int h_left = tmp.y;

	int tmp_h_left = h_left;

	#define s_do(x) single(h_max, f_top_##x, h_top_##x, e_left, h_left, h_topleft, score_##x);

	s_do(0); s_do(1); s_do(2); s_do(3); s_do(4); s_do(5); s_do(6); s_do(7);

	#undef s_do

	h_topleft = tmp_h_left;

	tmp.x = e_left;
	tmp.y = h_left;

	eh_tmp[offset] = tmp;
};

__device__ void load_letters(int k, int i, int step, uchar4* base_u, int& letter_0, int& letter_1, int& letter_2, int& letter_3, int& letter_4, int& letter_5, int& letter_6, int& letter_7){

	int offset = 2 * k * step + i;
	int offset2 = offset + step;

	uchar4 letter_container = base_u[offset];
	uchar4 letter_container2 = base_u[offset2];

	letter_0 = letter_container.x;
	letter_1 = letter_container.y;
	letter_2 = letter_container.z;
	letter_3 = letter_container.w;
	letter_4 = letter_container2.x;
	letter_5 = letter_container2.y;
	letter_6 = letter_container2.z;
	letter_7 = letter_container2.w;
};

__device__ void read_helper_data(int& offset, int step, int* base_helper, int& f_top_0, int& h_top_0, int& f_top_1, int& h_top_1, int& f_top_2, int& h_top_2, int& f_top_3, int& h_top_3, int& f_top_4, int& h_top_4, int& f_top_5, int& h_top_5, int& f_top_6, int& h_top_6, int& f_top_7, int& h_top_7){

	short4 f_data = (*((short4**)&(base_helper)))[offset];

	offset += step;
	short4 h_data = (*((short4**)&(base_helper)))[offset];

	offset += step;
	short4 f_data_2 = (*((short4**)&(base_helper)))[offset];

	offset += step;
	short4 h_data_2 = (*((short4**)&(base_helper)))[offset];
	offset += step;

	f_top_0 = f_data.x;
	h_top_0 = h_data.x;
	f_top_1 = f_data.y;
	h_top_1 = h_data.y;
	f_top_2 = f_data.z;
	h_top_2 = h_data.z;
	f_top_3 = f_data.w;
	h_top_3 = h_data.w;

	f_top_4 = f_data_2.x;
	h_top_4 = h_data_2.x;
	f_top_5 = f_data_2.y;
	h_top_5 = h_data_2.y;
	f_top_6 = f_data_2.z;
	h_top_6 = h_data_2.z;
	f_top_7 = f_data_2.w;
	h_top_7 = h_data_2.w;
};

__device__ void init_helper_data(int& offset, int step, int* base_helper, int& f_top_0, int& h_top_0, int& f_top_1, int& h_top_1, int& f_top_2, int& h_top_2, int& f_top_3, int& h_top_3, int& f_top_4, int& h_top_4, int& f_top_5, int& h_top_5, int& f_top_6, int& h_top_6, int& f_top_7, int& h_top_7){

	f_top_0 = h_top_0 = f_top_1 = h_top_1 = f_top_2 = h_top_2 = f_top_3 = h_top_3 = f_top_4 = h_top_4 = f_top_5 = h_top_5 = f_top_6 = h_top_6 = f_top_7 = h_top_7 = 0;
	offset += 4 * step;
};

__device__ void dump_helper_data(int& offset, int step, int* base_helper, int& f_top_0, int& h_top_0, int& f_top_1, int& h_top_1, int& f_top_2, int& h_top_2, int& f_top_3, int& h_top_3, int& f_top_4, int& h_top_4, int& f_top_5, int& h_top_5, int& f_top_6, int& h_top_6, int& f_top_7, int& h_top_7){

	short4 f_data, f_data_2, h_data, h_data_2;

	f_data.x = f_top_0;
	h_data.x = h_top_0;
	f_data.y = f_top_1;
	h_data.y = h_top_1;
	f_data.z = f_top_2;
	h_data.z = h_top_2;
	f_data.w = f_top_3;
	h_data.w = h_top_3;

	f_data_2.x = f_top_4;
	h_data_2.x = h_top_4;
	f_data_2.y = f_top_5;
	h_data_2.y = h_top_5;
	f_data_2.z = f_top_6;
	h_data_2.z = h_top_6;
	f_data_2.w = f_top_7;
	h_data_2.w = h_top_7;

	offset -= 4*step;

	(*((short4**)&(base_helper)))[offset] = f_data;

	offset += step;
	(*((short4**)&(base_helper)))[offset] = h_data;

	offset += step;
	(*((short4**)&(base_helper)))[offset] = f_data_2;

	offset += step;
	(*((short4**)&(base_helper)))[offset] = h_data_2;
	offset += step;
};

__device__ void do_strip(int z, short2* eh_tmp, int i, int block_len, int step, uchar4* base_u, int& h_max, int* base_helper, int score[shared_strip_width][alphabet_len]){

	for (int j = 0; j < shared_strip_width; j++){ 
		eh_tmp[BLOCK_SIZE * j] = make_short2(0, 0);
	};

	__syncthreads();

/*	if (threadIdx.x < shared_strip_width){
		for (int k = 0; k < alphabet_len; k++){
			score[threadIdx.x][k] = profile_matrix[z + threadIdx.x][k];
		};
	};*/

	if (threadIdx.x < alphabet_len){
		for (int k = 0; k < shared_strip_width; k++){
			score[k][threadIdx.x] = profile_matrix[z + k][threadIdx.x];
		};
	};


	__syncthreads(); 

	int h_toplast = 0;

	int off = i;

	for (int k = 0; k < block_len/8; k++){ // <- zwrocic na to uwage
		int letter_0, letter_1, letter_2, letter_3, letter_4, letter_5, letter_6, letter_7;
		load_letters(k, i, step, base_u, letter_0, letter_1, letter_2, letter_3, letter_4, letter_5, letter_6, letter_7);

		int f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7;

		read_helper_data(off, step, base_helper, f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7);

		int h_topleft = h_toplast;
		h_toplast = h_top_7;

		#define UNROLL(x) do_8step(eh_tmp, h_topleft, h_top_0, f_top_0, h_top_1, f_top_1, h_top_2, f_top_2, h_top_3, f_top_3, h_top_4, f_top_4, h_top_5, f_top_5, h_top_6, f_top_6, h_top_7, f_top_7, h_max, score[x][letter_0], score[x][letter_1], score[x][letter_2], score[x][letter_3], score[x][letter_4], score[x][letter_5], score[x][letter_6], score[x][letter_7], BLOCK_SIZE * x);

		#include __UNROLL_FILE__

		#undef UNROLL 

		dump_helper_data(off, step, base_helper, f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7);
	};
};

__device__ void do_init_strip(int z, short2* eh_tmp, int i, int block_len, int step, uchar4* base_u, int& h_max, int* base_helper, int score[shared_strip_width][alphabet_len]){

	for (int j = 0; j < shared_strip_width; j++){ 
		eh_tmp[BLOCK_SIZE * j] = make_short2(0, 0);
	};

	__syncthreads();

	/*if (threadIdx.x < shared_strip_width){
		for (int k = 0; k < alphabet_len; k++){
			score[threadIdx.x][k] = profile_matrix[z + threadIdx.x][k];
		};
	};*/

	if (threadIdx.x < alphabet_len){
		for (int k = 0; k < shared_strip_width; k++){
			score[k][threadIdx.x] = profile_matrix[z + k][threadIdx.x];
		};
	};

	__syncthreads(); 

	int h_toplast = 0;

	int off = i;

	for (int k = 0; k < block_len/8; k++){ // <- zwrocic na to uwage
		int letter_0, letter_1, letter_2, letter_3, letter_4, letter_5, letter_6, letter_7;
		load_letters(k, i, step, base_u, letter_0, letter_1, letter_2, letter_3, letter_4, letter_5, letter_6, letter_7);

		int f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7;

		init_helper_data(off, step, base_helper, f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7);

		int h_topleft = h_toplast;
		h_toplast = h_top_7;

		#define UNROLL(x) do_8step(eh_tmp, h_topleft, h_top_0, f_top_0, h_top_1, f_top_1, h_top_2, f_top_2, h_top_3, f_top_3, h_top_4, f_top_4, h_top_5, f_top_5, h_top_6, f_top_6, h_top_7, f_top_7, h_max, score[x][letter_0], score[x][letter_1], score[x][letter_2], score[x][letter_3], score[x][letter_4], score[x][letter_5], score[x][letter_6], score[x][letter_7], BLOCK_SIZE * x);

		#include __UNROLL_FILE__

		#undef UNROLL

		dump_helper_data(off, step, base_helper, f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7);
	};
};

template<int leftover>
__device__ void do_last_strip(int z, short2* eh_tmp, int i, int block_len, int step, uchar4* base_u, int& h_max, int* base_helper, int score[shared_strip_width][alphabet_len]){

	for (int j = 0; j < shared_strip_width; j++){ 
		eh_tmp[BLOCK_SIZE * j] = make_short2(0, 0);
	};

	__syncthreads();

/*	if (threadIdx.x < shared_strip_width){
		for (int k = 0; k < alphabet_len; k++){
			score[threadIdx.x][k] = profile_matrix[z + threadIdx.x][k];
		};
	};*/

	if (threadIdx.x < alphabet_len){
		for (int k = 0; k < shared_strip_width; k++){
			score[k][threadIdx.x] = profile_matrix[z + k][threadIdx.x];
		};
	};


	__syncthreads(); 

	int h_toplast = 0;

	int off = i;

	for (int k = 0; k < block_len/8; k++){ // <- zwrocic na to uwage
		int letter_0, letter_1, letter_2, letter_3, letter_4, letter_5, letter_6, letter_7;
		load_letters(k, i, step, base_u, letter_0, letter_1, letter_2, letter_3, letter_4, letter_5, letter_6, letter_7);

		int f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7;

		read_helper_data(off, step, base_helper, f_top_0, h_top_0, f_top_1, h_top_1, f_top_2, h_top_2, f_top_3, h_top_3, f_top_4, h_top_4, f_top_5, h_top_5, f_top_6, h_top_6, f_top_7, h_top_7);

		int h_topleft = h_toplast;
		h_toplast = h_top_7;

		#define DO_STEP(x) do_8step(eh_tmp, h_topleft, h_top_0, f_top_0, h_top_1, f_top_1, h_top_2, f_top_2, h_top_3, f_top_3, h_top_4, f_top_4, h_top_5, f_top_5, h_top_6, f_top_6, h_top_7, f_top_7, h_max, score[x][letter_0], score[x][letter_1], score[x][letter_2], score[x][letter_3], score[x][letter_4], score[x][letter_5], score[x][letter_6], score[x][letter_7], BLOCK_SIZE * x);

		for (int i = 0; i < leftover; i++){
			DO_STEP(i);
		};
		#undef DO_STEP 
	};
};

template<int leftover>
__global__ void psi_sw_block_3_template(int* base, int* base_helper, int* base_r, int block_len, int query_len){

	const int i = threadIdx.x + blockIdx.x * blockDim.x;
	const int step = blockDim.x * gridDim.x;

	const int buf_len = shared_strip_width * BLOCK_SIZE;

	__shared__ short2 eh_tmp_base[buf_len];
	short2* eh_tmp = eh_tmp_base + threadIdx.x;
	uchar4* base_u = (uchar4*)base;

	__shared__ int score[shared_strip_width][alphabet_len];

	int h_max = 0;
	int z = 0;

	do_init_strip(z, eh_tmp, i, block_len, step, base_u, h_max, base_helper, score);

	z += shared_strip_width;

	while (z < query_len){

		do_strip(z, eh_tmp, i, block_len, step, base_u, h_max, base_helper, score);
				
		z += shared_strip_width;
	};

	do_last_strip<leftover>(z, eh_tmp, i, block_len, step, base_u, h_max, base_helper, score);

	base_r[i] = h_max;
};

};

#include <stdio.h>

/** 
  Computation

  One should init gap_open, gap_extend and fill in profile_matrix matrix with proper profile matrix (see below).
*/

/**
  Database block format.

  This kernel processes BLOCK * BLOCK_SIZE sequences at once. All sequences are block_len
  (see run_psi_sw_block_3_new arguments) long or are padded with neutral symbol. Amino acids
  are represented by codes that are 4 bytes long (int). We use values from 0 to alphabet_len only.

  Profile matrix is constructed by copying columns with scores for amino acids in query from ordinary
  score matrix and inserting neutral score (which is 0) at position we choose. When generating
  data for GPU to process we need to pad sequences that are too short with number of this position.
*/


/** 
  @param BLOCKS number of blocks to run, for optimal execution performance it should be one per multiprocessor
  @param base pointer to gpu memory with database block to process. BLOCK_SIZE * BLOCKS * block_len * sizeof(int)
  @param base_helper pointer to gpu memory to use for auxiliary data
  @param base_r pointer to gpu memory for result storing
  @param block_len base length
  @param query_len query length
  @param stream ignored parameter
*/

extern "C" void run_psi_sw_block_3_new(int BLOCKS, void* base, void* base_helper, void* base_r, int block_len, int query_len, cudaStream_t stream){

	int q_org = query_len;

	if (query_len % new_hope_3::shared_strip_width != 0){
		query_len -= query_len % new_hope_3::shared_strip_width;
	};

	#define UNROLL(x) case x: { new_hope_3::psi_sw_block_3_template<x><<<BLOCKS, new_hope_3::BLOCK_SIZE, 0>>>((int*)base, (int*)base_helper, (int*)base_r, block_len, query_len); } break;

	int switch_value = q_org % new_hope_3::shared_strip_width;

	if (q_org < new_hope_3::shared_strip_width){
		switch_value = 0;
	};

	switch (switch_value){

		#include __UNROLL_FILE__
	};

	#undef UNROLL
};

extern "C" int strip_width(){
	return new_hope_3::shared_strip_width;
}

extern "C" int block_size(){
	return new_hope_3::BLOCK_SIZE;
}

#undef __UNROLL_FILE__

#endif
